Glass transition and layering effects in confined water: 
a computer simulation study 
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Single particle dynamics of water confined in a nanopore is studied through Computer Molecular 
Dynamics. The pore is modeled to represent the average properties of a pore of Vycor glass. 
Dynamics is analyzed at different hydration levels and upon supercooling. At all hydration levels 
and all temperatures investigated a layering effect is observed due to the strong hydrophilicity of the 
substrate. The time density correlators show, already at ambient temperature, strong deviations 
from the Debye and the stretched exponential behavior. Both on decreasing hydration level and 
upon supercooling we find features that can be related to the cage effect typical of a supercooled 
liquid undergoing a kinetic glass transition. Nonetheless the behavior predicted by Mode Coupling 
Theory can be observed only by carrying out a proper shell analysis of the density correlators. 
Water molecules within the first two layers from the substrate are in a glassy state already at 
ambient temperature (bound water). The remaining subset of molecules (free water) undergoes a 
kinetic glass transition; the relaxation of the density correlators agree with the main predictions of 
the theory. From our data we can predict the temperature of structural arrest of free water. 
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I. INTRODUCTION 



The modification of the dynamical properties of liq- 
uids in confined environment relative to the bulk is a 
field of rapidly growing interest because of the close con- 
nection with relevant technological and biophysical prob- 
lems. The specific differences in behavior are, among 
others, due to different interactions between liquid and 
substrate, the size of the confining region, and the size 
of the particles composing the liquid. Nonetheless some 
underlying common features can be sorted out of thej-ex- 
tremely rich and diversified phenomenology available .u In 
fact there are two main causes for the change of dynam- 
ics of a confined liquid with respect to its bulk phase: 
the first is the influence of the interactions of the liquid 
with a rough surface, which is expected to slow down 
the dynamics; the second is the confinement effect that 
can lead to an increase in the free volume of a molecule, 
which then results in accelerated dynamics together with 
a decrease of the glass transition temperature relative to 
the one of the bulk phase. The interplay between these 
two effects depends crucially onJjhe particle density and 
the size of the confining system.cTEl 

In particular the dynamical properties of water in re- 
stricted geometries and at interfaces have recently been 
studied intensely because of the important effects in sys- 
tems of interest to biology, chemistry, and geophysics, 



whose behavior depends on how the pore size and struc- 
ture influence the diffusion of water. Those properties are 
particularly relevant in understanding phenomena like 
the mobility of water in biological channelscrQ or the dy- 
namics of hydrated proteins Jj 

It is well known that liquid water shows a very peculiar 
behavior in the supercooled phase. The study of water 
approaching a glass phase is still a challenging problem 
since below 235 K one enters the so called no man's land 

where nucleation processes take place and drive the 
liquid to the solid crystalline phase, preventing the obser- 
vation of the glass transition.E30 It is unclear until now 
how confinement could change this scenario. It would be 
particularly important to understand whether the glass 
transition temperature could be experimentally accessi- 
ble for confined water. In this respect the modifications 
induced by the confinement on the dynamics of water on 
supercooling are of extreme interest. 

Computer simulation is a very suitable tool for ex- 
ploring the liquid in the range of the supercooled regime 
without the limitations of the nucleation process which 
takes place in the real experiment. Dynamical proper- 
ties of different types of liquids in confined geometries 
have|-|bae&, studied in recent years by computer simula- 
tion.BOiiJ For water confined in micropores there are a 
number-of computer simulation studies on the mobility of 
water .Bl3't£rE3 It is still difficult, however, to find general 
trends; systematic studies of the dynamics of confined 
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water have not been attempted until now. 

For confined and interfacial water inelastic neutron 
scattering and NMR spectroscopy found a slowing down 
of the dynamics relative to the bulk phaseEOiSil For wa- 
ter in contact with proteins there are also signatures of 
a transition of adsorbed water-to a glassy state, which is 
driven by the protein surface!!! Moreover, recent simula- 
tion and experimental studies showed a t ypic al spectral 
glassy anomaly, the so called boson peak.EaE3 

In the many experimental studies of confined water 
that has*! been performed water in Vycor is of particular 
interestLD, since Vycor is a porous silica glass charac- 
terized by a quite sharp distribution of pore sizes with 
an average diameter of 40 A. The pore size does not de- 
pend on the hydration level and the surface of the pore is 
strongly hydrophilic. Moreover, the water-in- Vycor sys- 
tem can be considered as a prototype for more complex 
situations of interfacial water .EJ'Ej 

In this paper we present the results of a Molecular Dy- 
namics study of the single particle dynamics of water con- 
fined in a silica cavity modeled to rep re s en t the average 
properties of the pores of Vycor glass.L§L3 We will con- 
centrate on the dynamical behavior of the confined water 
at half hydration on supercooling. In previous studie^-the 
static properties of this system were investigated.c3E3 

In the second section we recall the general theoretical 
background which concerns the modern interpretation of 
the glass transition and we will concentrate on a.short de- 
scription of the Mode Coupling Theory (MCT). El In sec- 
tion III we explain the technical details of our simulation 
of Molecular Dynamics (MD). In section IV we discuss 
the important layering effect that we observe in the con- 
fined water. This effect forms the basis for a layer anal- 
ysis of the dynamieal results that we introduced for the 
density correlator B3 In Section V this analysis is carried 
out for the system at ambient conditions. By means of 
this layer analysis we can find agreement with the MCT 
predictions and estimate the values of some relevant pa- 
rameters of the theory as exposed in Section VI where 
the properties of the system upon supercooling are ana- 
lyzed. Section VII is devoted to the conclusions of our 
work. 



A scenario emerges that can be pictured by considering 
the free energy landscape or alternatively the potential 
energy landscapes In the high temperature regime the 
free energy functional has only one relevant minimum and 
correspondingly there are several directions in the con- 
figurational space with negative eigenvalues. The liquid 
behaves "normally", and diffusion is stochastic, i.e., there 
is only one relevant timescale and the relaxation process 
is of Debye type. As one starts supercooling (i.e. cool- 
ing the system below the melting temperature T m ) the 
system approaches an important crossover temperature 
Tc referred to. as the temperature of kinetic glass tran- 
sition. MCTc3 in its idealized version is able to describe 
the dynamics of the liquid in great detail with precise 
predictions on the behavior and the analytical shape of 
the density correlators (see next subsection) for systems 
sufficiently close to, yet above-jTc- The region above Tc 
is only landscape influenced. l3 At T = Tc the nup 
of directions leading to different basins goes to zero.l 
In the region below Tc there is an exponentially large 
number of minima in the free energy separated by bar- 
riers that are high compared to thermal energy. In this 
respect Tc represents a cross-over temperature from a 
liquid-like regime to a solid-like regime where only hop- 
ping processes can restore ergodicity. Below Tc the su- 
percooled liquid is in fact frequently trapped in one of the 
local free energy minima and only activated processes 
take the system frorn-pne minimum to another one in 
the energy landscape^ This region is strongly landscape 
dominated. L3 

In the "many valley picture" the liquid-glass transition 
occurs at the Kauzmann temperature Tk, where the con- 
figurational entropy, which measures the logarithm of the 
number of accessible minima, vanishes .and, .correspond- 
ingly, the free energy barriers diverge .EErtiloa The tran- 
sition at Tk can be considered as an ideal glass transition 
which can take place only at an infinitely slow cooling 
rate, where it would be signaled by the divergence of the 
viscosity. This ideal second order transition is related to 
the singularities which are found at finite cooling rate in 
experiments at the conventional glass transition temper- 
ature, T g , where Tk <T g < Te- 



ll. THE PHYSICS OF THE GLASS TRANSITION. 

A. The glass transition scenario 

Glasses are a very popular subject of the contemporary 
science literature. But in spite of the huge experimental, 
computational and theoretical efforts made, a unified pic- 
ture able to account for the behavior of the liquid from 
its normal state throughout supercooling and vitrifica- 
tion is still lacking. Theoretical physicists are develop- 
ing a first principle theory for glassesE3Ea, based pa the 
phenomenological ideas introduced by Kauzmanetp and 
developed later by Adam, Gibbs and Di MarzioS 



B. Predictions of the Mode Coupling Theory 

In this paper we will focus on the relatively high tem- 
perature region where dynamics can be studied by MD, 
i.e. the supercooled region where T is above and ap- 
proaches Tc- Here MCT io_its idealized version works 
very well for many systems £11 

MCT is able to describe the dynamics of a liquid when 
the single entity, molecule or atom, is trapped by the 
transient cage formed by its nearest neighbors. This 
transient caging is responsible for the stretching of the 
relaxation laws and the separation of time scales. In its 
idealized version MCT does not take into account the 
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hopping processes. This version predicts a transition to 
a non-ergodic system at Tc, when all the cages are frozen. 
Nonetheless, hopping processes are not relevant for most 
liquids above Tc- An asymptotic expansion in the re- 
gion near the ideal singularity (at Tc) yields predictions 
of the functional form of the time correlators and the 
critical exponents. These predictions have been tested 
both by experiments and MD simulations for many glass 
formers. MCT thus renders possible, among others, an 
estimate of the behavior o f the time correlators close to 
the crossover temperature. Eijcil 

MCT predicts a two-step relaxation for the dynam- 
ics, which is the signature of the cage effect. Well above 
Tc, after an initial ballistic regime, the particle enters 
the stochastic diffusion regime, and its mean square dis- 
placement (MSD) increases linearly with time. In the su- 
percooled region above Tc, after the short time ballistic 
regime, the particle is trapped by the barrier created by 
the nearest-neighbor cage; Brownian diffusion is restored 
only when this cage relaxes. This behavior is reflected 
in the self part of the density autocorrelation function, 
and in any other correlator which has a nonzero overlap 
with the density. In particular, when approaching Tc 
from the liquid side, the Fourier transform of the den- 
sity correlator, the self intermediate scattering function 
(ISF), Fs(Q,t), has a two step relaxation behavior with 
a fast and a slow decay. After the fast decay it enters a 
plateau region, corresponding to the rattling of the par- 
ticle in the nearest-neighbor cage, which is called the j3- 
relaxation region. After the time interval of the plateau, 
which becomes longer when approaching Tc, the func- 
tion Fs(Q,t) decays to zero. This long time relaxation 
is called a-relaxation. In the a-relaxation region it has 
been shown that the relaxation process is well described 
by a stretched exponential 
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where 77 is the long relaxation time and /3 is called the 
Kohlrausch exponent. MCT predicts that, when Tc is 
approached, the a-relaxation takes place on increasingly 
longer time scales, so that the relaxation time 77 diverges 
with a power law 



(T - T c ) 



(2) 



As a consequence, the diffusion coefficient D, which is 
predicted to be proportional to t^ 1 , goes to zero at Tc 
with the power law D ~ (T- Tc) 1 ■ 

In the idealized version of MCT the system becomes 
non-ergodic at Tc, defined as the temperature of struc- 
tural arrest; below Tc the cages are frozen and only hop=. 
ping processes can restore ergodicity. Extended MCTc3 
takes into account this effect. 

Finally, we recall that, since the MCT description of 
the kinetic glass transition is based on the cage effect, 
the relevant length scales are of the order of the nearest- 
neighbor distances. Consequently the dynamical quanti- 
ties in Q space display this effect most clearly for values 



of Q close to the maximum of the static structure factor 
S(Q). 



C. Bulk supercooled water and MCT 

In recent years it was discovered that spc/eEI super- 

Boled water has a temperature of structural arrest Tc 
coinciding with the_sp called singular temperature T s 
of Speedy and Angell.c3 As stated above, MCT predicts 
that close to Tc the liquid dynamics is dominated by the 
"cage effect" . Water does not behave like normal glass- 
forming fluids in this regard, since the cage effect is not a 
consequence of an increase of density upon supercooling 
but rather seems to be determined by the increase of the 
hydrogen bond stiffness, which makes the cage more rigid 
as the temperature is lowered below room temperature. 
The long time behavior of the single particle dynamics 
is well described in terms of the MCT and the depen- 
dence of 77 (and D) on temperature are found to agree 
with the power law, Eq.(^J). Successive simulations over 
a wide range of-pressures and temperatures^ and the- 
oretical studiesc3 fully confirmed the MCT behavior of 
this potential. 

TABLE I. Parameters of the SPC/E model and of the sub- 
strate-water interaction potential. 





Site 


a 

( A ) 


(K) 


q 

|e| 


water a 





3.154 


78.0 


-1.04 




H 


0.0 


0.0 


0.52 




H 


0.0 


0.0 


0.52 


silica b 


Si 


0.0 


0.0 


1.283 




BO 


2.70 


230.0 


-0.629 




NBO 


3.00 


230.0 


-0.533 




AH 


0.0 


0.0 


0.206 


a SPC/E model.@ 
b Values from Ref. 
ing oxygens; AH: 


[si|; BO: bridging oxygens; NBO: non 
protons on the substrate surface. 


bridg- 



TABLE II. Hydration levels of the pore. Nw is the num- 
ber of water molecules and pw the corresponding global den- 
sity. The hydration level is based on estimated value for full 
hydration Nty ~ 2600 molecules (see text). 



N 



% hydration 



pw {g/cm 



500 
1000 
1500 
2000 
2600 



19% 
38% 
56% 
75% 
98% 



0.1687 
0.3373 
0.4971 
0.6658 
0.8788 
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III. MOLECULAR DYNAMICS OF WATER IN 
VYCOR: TECHNICAL DETAILS 



Vycor is a porous silica glass obtained by spinodal de- 
composition of a glass- forming melt of SiC>2 and B2O3. 
The B2C>3-rich phase is leached out, leaving a SiC>2 
porous glass. The main differences of Vycor glass with 
respect to other analogous porous media are that it has 
a well characterized structure with a quite sharp distri- 
bution of pore sizes and a strong hydropylic surface. The 
void fraction, 28%, is an interconnected network of pores 
of diameter ~ 40 A. Vycor has a strong capability to 
absorb water, since its equilibrium hydration at ambi- 
ent conditions is 25% of its dry weight. Furthermore, its 
pores do not change size when filled with water. For these 
reasons it can be considered as a good candidate to study 
the general behavior of water in hydcaphilic nanopores as 
a function of the level of hydration.o 

During experimental sample preparation, after the 
dessication process and before introducing water into the 
pores, those oxygen atoms on the internal surfaces whose 
valences are not saturated (called non-bridging oxygens) 
are saturated with hydrogen atoms. 

In our simulation we build up a cubic cell of silica glass 
by the usual procedure of melting a /3-cristobalite crys- 
talline structure at 6000 K and quenching to room, tem- 
perature. As described in details in previous workOEZI we 
get a cube of length L = 71.29 A. Inside the cube we carve 
a cylindrical cavity of 40 A diameter by eliminating all 
the atoms lying within a distance R — 20 A from the axis 
of the cylinder taken as the z-axis. Then we exclude all 
the silicons with less than four oxygen neighbors. On the 
surface of the pore we distinguish the oxygen atoms which 
are bonded to two silicons (so-called bridging oxygens 
(BO)) from those which are bonded to only one silicon 
atom called non-bridging oxygens (NBO). The NBOs are 
saturated with hydrogen atoms. The final block of ma- 
terial contains 6155 silicon atoms, 12478 BOs, 227 NBOs 
saturated by 227 acidic hydrogens (AHs). The surface 
density of AH results to be 2.5 nm~ 2 , in good agreement 
with the experimentally determined value of 2.3 nm~ 2 .n 

Water molecules described by— .the Simple Point 
Charge/Extended (SPC/E) modeled are introduced in 
the cavity. The SPC/E potential is one of the most 
frequently used interaction site models for water. Both 
density and diffusion coefficient show a remarkably good 
agreement with experiment at ambient conditioned It 
has to be stressed that SPC/E exhibits a minimum in 
the pressure P versus T at about 240 K, while the ex- 
perimental temperature of maximum density is located 
at T ~ 277 K. We therefore expect the phase diagram 
of SPC/E water to be shifted downwards in temperature 
relative to the experimental one.c3l£3 

The water sites interact with the atoms of the xigid ma- 
trix by means of an empirical potential modeEj'ELj, where 
different fractional charges are assigned to the sites of the 
silica glass. In addition, the oxygen sites of water interact 



with BOs and NBOs of the substrate via Lennard- Jones 
potentials. All parameters are collected in Table |. 

The molecular dynamics is performed in the micro- 
canonical ensemble with a timestep of 2.5 fs. Each run 
was equilibrated via a coupling to a temperature "reser- 
voir" by using the Berendsen method of velocity rescal- 
ing. Data during the initial equilibration period were 
discarded, until monitored quantities like the internal en- 
ergy showed no trends with time. 

Since the simulations of such a large system are rather 
time consuming, we used the shifted force method with 
a cutoff at 9 A for all interactions. We checked that 
the use of a larger cutoff or Ewald summations does not 
change the tread of the obtained results, as discussed in 
previous work.EH Since we need to store long trajectories 
for a large number of molecules, we minimized the disk 
space by saving them at a logarithmic time step. With 
this choice the size of trajectories for 2000 molecules last- 
ing 1 ns is roughly 1.5 GB. For the lower temperatures 
several runs of 1 ns each were performed. 

In the following we will present data for five hydration 
levels of the pore at ambient temperature, as described 
in table O. For the case of roughly half hydration, which 
corresponds to Nyy — 1500 water molecules, we stud- 
ied five temperatures, namely T = 298, 270, 240, 220 and 
210 K. 
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FIG. 1. Snapshots of equilibrium configurations of confined 
water at different hydration levels (see Table II). Only the 
oxygen atoms of water are shown projected on the xy plane 
perpendicular to the axis of the confining cylinder. 



IV. DENSITY PROFILES OF CONFINED WATER 

In the experiments on water confined in Vycor the full 
hydration is obtained by immersion of the dry sample in 
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water. Vycor glass adsorbs water up to 25% of its dry 
weight, so the full hydration in the experiment is defined 
as hf ~ 0.25 g of water/g of Vycor. A partial hydration 
is experimentally obtained by exposing the full hydrated 
sample to P2O5 for a number of hours in order to re- 
duce the water content and reach the desired ratio h — 
mass of water/mass of Vycor. c3 

The full hydration in the_experiment (h = hf) corre- 
sponds to a density of waterEj pf = 0.8877 gcm~ 3 , which 
is 11% lower respect to bulk water. In computer simula- 
tion we vary the hydration level by changing the number 
of molecules inserted in the single cylindrical cavity of our 
simulation cell. We assume that the full hydration cor- 
responds to the experimental density pf which would be 
given in our simulation by a number of water molecules 
Nw ~ 2600. Levels of hydration are given relative to 
this value. 



(~ 20%), N w = 1000 (~ 40%), N w = 1500 (~ 60%), 
N w = 2000 (~ 75%) and N w = 2600 (~ 100%), as 
in Table ||. We observe that the pore is strongly hy- 
drophilic, since molecules bind strongly to the surface 
even at low levels of hydration. The lowest hydration 
level is lower than the .monolayer coverage estimated 
to be 25% of hydration.^ In the inner surface of the 
pore single water molecules or small patches of water 
molecules are found. The radial density profiles normal- 
ized to the bulk water density are shown in Fig. || as a 
function of R = y/ ' x 2 + y 2 . The formation of well de- 
fined layers of water molecules close to the substrate is 
evident in the region between 15 A and 20 A, where the 
pore surface is located. For the highest level of hydration 
water fills the pore completely. In the pore center a den- 
sity close to the experimental density at full hydration is 
reached. 
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FIG. 2. Density pr ofiles of confined water along the cylin- 
drical radius R = v/x 2 + y 2 normalized to the bulk water 
density at ambient temperature. Top: density profiles at 
T — 298 for decreasing levels of hydration (from Nw = 2600 
to Nw = 500 from top to bottom). A small fraction of wa- 
ter molecules is trapped inside the silica glass, which leads to 
the density contribution for R > 20 A, where the confining 
surface is located. Bottom: density profiles for Nw = 1500 
at ambient temperature, T = 298 K, (continuous line) and in 
the supercooled regime, T = 210 K, (long dashed line). 

In Fig. [l] we show snapshots at ambient temperature 
for the five different hydration levels studied, Nw = 500 
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FIG. 3. Mean square deviation (MSD) of the oxygen atoms 
< r 2 {t) >=< \r{t) -r(0)| 2 > for diffusion in three dimen- 
sions at ambient temperature as a function of time, t, for 
different levels of hydration from top {Nw = 2600) to bottom 
{Nw = 500). In the inset the diffusion coefficients extracted 
from the slopes of the MSD curves are displayed as a function 
of the hydration level (filled circles). For comparison the bulk 
value for SPC/E water at ambient conditions is also shown 
(filled diamond). 



Layering effects of confined water have been observed 
in almost all computer simulations performed in differ- 
ent geometry and with different water-substrate inter- 
actionJla In our case they are representative of the high 
hydrophilicity of the pore. In particular, we note that, as 
the hydration level is increased, a double layer structure 
with density oscillations is formed. Another consequence 
of the hydrophilic interaction is the distortion of the hy- 
drogen bo nd n etwork of water in the layers close to the 
substrate.EScZI Since the inner surface is corrugated, we 
do observe that the density is different from zero also 
right inside the pore surface, where water molecules are 
trapped in small pockets inside the substrate. These 
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molecules do not contribute to the diffusion dynamics 
of water. A decrease of temperature has no large effect 
on the density profiles, as can be seen in the lower part 
of Fig. H where the case Nw = 1500 is shown at ambient 
temperature (T = 297 K) and in the supercooled regime 
(T = 210 K). We note that the peak heights of the two 
layers become more pronounced as the temperature is 
lowered. 



V. SLOW DYNAMICS OF CONFINED WATER: 
THE ROLE OF HYDRATION AT AMBIENT 
TEMPERATURE 

We will consider now the single particle translational 
dynamics of the oxygens of water as a function of hy- 
dration level at ambient temperature. In Fig. we dis- 
play the mean square displacement (MSD) < r^t) >=< 
|r(t) — r(0)| 2 > in three dimensions < r 2 >—< x 2 > + < 
y 2 > + < z 2 > as a function of time for the different hy- 
dration levels. The plot is on a log-log scale in order to 
make more evident the flattening of the curves at lower 
hydrations. After the initial ballistic regime, where the 
MSD increases proportional to t 2 (up to t » 0.4 ps) the 
curves begin to flatten at intermediate times. This be- 
havior becomes more evident as the hydration level is 
decreased. While the onset of the cage effect is almost 
independent of the hydration, the relaxation time of the 
cage depends strongly on the hydration level. The lower- 
ing of the hydration affects the mobility of water already 
at room temperature by shifting the onseL-pf the diffu- 
sive behavior (MSD cx t) to larger times&Eil In the inset 
the diffusion coefficient, extracted from the slope of the 
MSD in the diffusive regime, is plotted as a function of 
the number of particles inside the pore. The SPC/E bulk 
water value is also displayed. We note substantial de- 
crease of the mobility of confined water as the hydration 
level is lowered. 

As can be expected, the different regimes in the diffu- 
sion are also reflected in the behavior of the single particle 
intermediate scattering function (ISF), Fs(Q,t). This is 
the Fourier transform of the Van Hove self-correlation 
function 

G S (r,t) = ^(j2S[r + r i (0)-r i (t)}\ (3) 

Gs(r, t)dr is proportional to the probability of finding a 
particle at distance r after a time t if the same particle 
was in the origin r = at the initial time t = 0. The 
incoherent, or self, ISF can be written as: 

Fs(Q,t) = (^e^ r ^- T ^ (4) 

In Fig. H we show the function Fs(Q,t) for the different 
hydrations calculated at the maximum of the oxygen- 
oxygen structure factor Q max — 2.25 A -1 , where the 



cage effect is expected to be strongest (see discussion in 
Sec. II B). A shoulder appears in the correlators upon 
decreasing the hydration, and the long time relaxation 
behavior deviates strongly from the exponential form. It 
is evident that the interaction with the hydrophilic sub- 
strate plays a crucial role in the cage effect, since, at 
the lower levels of hydration, the fraction of molecules in 
contact with the surface is larger than at higher levels of 
hydration. 
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FIG. 4. Self intermediate scattering function (ISF) of the 
oxygen atoms at ambient temperature for Q — 2.25 A -1 , the 
position of the maximum of the oxygen structure factor, for 
different levels of hydration from Nw = 500 (top curve) to 
Nw = 2600 (bottom curve). Q is averaged over all three 
dimensions. 

Apparently, lowering the hydration plays a role simi- 
lar to supercooling the bulk liquid, but, in spite of some 
similarity with bulk supercooled water, here the late part 
of the ISFs cannot be fitted by the stretched exponential 
of Eq. (|), which is predicted by MCT. 

Thus, at first glance our study of confined water ap- 
pears to indicate that the water mobility is reduced rel- 
ative to the bulk phase and that, as a function of hydra- 
tion, a diversification of relaxation times develops anal- 
ogous to supercooled systems undergoing a kinetic glass 
transition. Nonetheless the deviation for the late part of 
the correlator from the analytical predictions of MCT has 
to be explained for a full understanding of the dynamics. 

We must take into account in our analysis that the re- 
laxation of the cage around the molecules close to the 
substrate is determined by the strong hydrophilic effect 
which leads to the formation of the double layer struc- 
ture of water close to the surface (see Fig. ||). The way 
in which the water molecules arrange themselves close to 
the substrate as well as the extension of the substrate 
perturbation affect the dynamical properties. 
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A more detailed examination of the contributions to 
the single particle dynamics must therefore be performed 
by taking into account the possibility that the water 
molecules close to the Vycor surface could be much less 
mobile than the ones in the center of the pore, since the, 
former ones are strongly H-bonded to the substrateBS'EZl 

For this reason we performed a layer analysis separat- 
ing the contribution coming from the first two layers (the 
outer shells or outer layers) close to the substrate from 
the contribution coming from the remaining water in the 
inner shells (or inner layers) O The definition of these 
two subsets is inspired by the density profile (see fig.||). 
The outer shells are defined as 15 < R < 20 A while the 
inner shells are the remaining < R < 15 A. 
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FIG. 5. Layer analysis of the MSD of the oxygen atoms 
at ambient temperature. The contributions of the water 
molecules in the inner layers (free water) and in the outer 
layers (bound water) are presented at ambient temperature 
for Nw = 2600 (upper panel) and Nw = 1500 (lower panel). 
The continuous lines represent diffusion along the z direction, 
while the dashed lines are for the diffusion along the xy di- 
rection. 

In Fig. | we show the results of this analysis for the 
MSD at two different hydrations, roughly full hydration 
and half hydration. Since the fluid is confined in the 
x and y directions we calculate separately the quantity 
along the z direction and along the radial (or xy) di- 
rection R. Both quantities are multiplied by the proper 



factor in order to adjust them to the three dimensional 
scale for a direct comparison. 

It is clearly seen that the diffusion of the molecules 
in the outer shells is much slower relative to the ones 
in the inner shells. We also observe that the change of 
hydration does not have a strong influence on the outer 
shells. We find that the diffusion is slower along xy but 
the shape of the MSD does not change relative to the 
one in the z direction. In fact, the outer water molecules 
do not reach the diffusive regime during the observation 
time of Fig. 5 for both directions. 
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FIG. 6. Layer analysis of the self intermediate scattering 
function (ISF) at ambient temperature for Nw = 2600 at the 
peak of the structure factor. The ISF is shown along the z 
direction (top) and along the xy direction (bottom). In both 
directions the contributions from free and bound water are 
shown separately in addition to the total ISF. The free water 
contribution is fitted by Eq. (^) (long dashed line). 

The results of the layer analysis for the ISF at full hy- 
dration are shown in Fig. [(] along the z direction (top) 
and along the xy direction (bottom) . Water in the outer 
shells appears to be already in a glasslike state, since its 
ISF decays to zero over a much longer time scale. The 
behavior of the ISF from the inner layer contribution is 
different. We observe that the long time relaxation is 
characterized by a stretched exponential form, Eq. (Q). 
The onset of the a relaxation is preceded by the fast re- 
laxation, which is usually observed at short time both ex- 
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perimentally and in computer simulations of supercooled 
liquids. By taking into account the fast relaxation with 
an appropriate normalizing function in both directions 
the entire curve can be described very well by two relax- 
ation processes with different time scales 

F S (Q, t) = [1 - MQ)\ e- {t/Ts)2 + A(Q)e-<-W (5) 

where the stretched exponential contains the long relax- 
ation time t\ and the Kohlrausch exponent already 
introduced in Sec. II B. In this equation the normalizing 
function A(Q) is the Debye- Waller factor, which is also 
termed the Lamb-Mossbauer factor for the single particle 
motion. A(Q) — er a ® / 3 accounts for the cage effect, 
with a the effective cage radius. The short time function 
is written in terms of an exponential containing the short 
relaxation time t s . The values extracted from the fit are 
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FIG. 7. Like Fig. | but for N w = 1500. 

In Fig. the shell analysis for the ISF at Q m ax = 2.25 
A is shown at half hydration, (Nw = 1500). The shell 
analysis works well also for this level of hydration, and 
the curves obtained from the inner layer contribution can 
be fitted by the functional form Eq. d^). The resulting 
parameters are also reported in table III. 

The t parameters reported in table are similar to 
those of bulk water but (3 is much smallerES relative to the 



bulk where (3 = l£3 From the values of A(Q) the cage 
radius can be extracted. We obtain a ~ 0.5 A, again 
similar to the radius obtained for bulk water. 

We note in passing that if only one layer is considered 
for the bound water, namely 17 < R < 20 A, then the 
fit to the stretched exponential, of the remaining water 
is not completely satisfactory.E3 This confirms that the 
separation into the two subsets of water molecules chosen 
above is an appropriate representation of the two clearly 
distinct dynamical regimes coexisting in the liquid. It 
is unclear to what extent these effects would be visible 
in the experiments, where usually only the total ISF is 
measured. There are, however, experimental indications 
that ±»p,different subsets exist in H-bonded confined liq- 
uids .E3'E3 

We note that water shows an analogous behavior 
to what was found in mixtures, for example, with n- 
butoxyethanol. In this case dielectric relaxation studies 
show the existence of two kinds of water coexisting in tha 
mixture: the hydration water and the "bulk" water.EZI 
Similar behavior has also been reported in a theoreti- 
cal study of water in contact with a formic acid dimer, 
where significantly lower values of D have been observed 
in the first and second solvation region up to 4.5 A.E3 
Moreover, the dynamics of a complex liquid, o-terphenyl 
(OTP), confined in nanometer scale pores can also be 
interpreted in terms of a two layer model. Calorimet- 
ric measurements show in fact the existence of-,two glass 
transition temperatures for, this latter system.el _ _ 

Thermometric studiesEJC3, NMR spectroscopi£a~o, 
neutron diffractionEjS~E3, and X-ray diffractionEll can 
also be interpreted in terms of two types of water which 
are present in confining pores: free water ; which is in 
the middle of the pore, and bound water, which resides 
close to the surface. Free water is observed to freeze 
abruptly in the cubic ice structure with an hysteresis 
in the melting/freezing transition which decreases at de- 
creasing pore size. Bound water freezes gradually with- 
out any hysteris effect but it does not make any transition 
to an ice phase and it is called sometimes nonfreezable 
water.E3 Although all these experimental findings refer 
to static properties, we can reasonably identify the water 
molecules in the outer layers as bound and not freezable 
water while free water corresponds to the molecules in 
the inner layers of the simulated pore. 

TABLE III. Fit parameters of the lineshape, Eq. (|H|), to 
the ISF data at ambient temperature. 
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VI. CONFINED WATER AT SUPERCOOLED 
TEMPERATURES 



For the discussion of supercooled confined water we 
consider the case close to half hydration, Nw = 1500, 
which can be regarded as a reasonable compromise be- 
tween computational effort and performing the layer 
analysis with sufficient statistical accuracy. As stated 
above, the temperature effect on the structure of the dou- 
ble layer of molecules close to the substrate is small (see 
Fig. ||). Therefore, the layer analysis is done using the 
same separation, JZ. = 15 A, as for the room temper- 
ature simulations £3 The MSD along the z direction is 
shown for decreasing temperatures for the inner shells 
(free water) in Fig. || and the outer shells (bound wa- 
ter) in Fig. |^. The clear distinction between free and 
bound water at room temperature is maintained at all 
lower temperatures. It is clear from Fig. ^ that at lower 
temperatures the cage effect persists for longer times. 
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FIG. 9. MSD of bound water at half hydration along the z 
direction for 5 temperatures (298, 270, 240, 220, and 210 K, 
from top to bottom). 
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FIG. 8. MSD of free water at half hydration along the z 
direction for 5 temperatures (298, 270, 240, 220, and 210 K, 
from top to bottom). 



A. The ISF at Q max versus T 



Fig. [lO] shows, as an example, the ISF for motion along 
the z direction at T = 240 K and Q = Q max = 2.25 A. 
The inner layer contribution can be fitted very well by 
Eq. (||); we obtain j3 — 0.62, r ; = 11 ps, t s = 0.16 ps. 
At this temperature a bump appears in the ISF of the 
free water at around 0.7 ps. Signatures of this feature 
are also visible in the MSD of Fig. || and ^ where oscil- 
lations in the plateau region are easily discernible. They 
are similar to features found in other simulation studies; 
for typically strong glass formers th | ey ,| Ca p-|be attributed 
to the so called Boson peak (BP).E3He§ti2l 
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FIG. 10. Layer analysis of the ISF at T = 240 K for 
N w = 1500 at Q = 2.25 A" 1 along the z direction. The 
contributions from the free (inner layers) and bound water 
(outer layers) are shown together with the total ISF. The free 
water curve is fitted by Eq. (||) (long dashed line). The bump 
in the free water ISF is evidence of the presence of a Boson 
Peak. 

We repeated the layer analysis of the ISF at Q max for 
all different temperatures studied. In Fig. |ll] we show 
the results along the z direction for the inner layers as 
a function of temperature at Q = Q ma x- We wish to 
point out that there are no qualitative differences in the 
behavior along the xy direction, where, in particular, the 
same BP feature is present. We observe that, as the tem- 
perature is decreased, the shoulder of the slow relaxation 
becomes more and more pronounced. 







The fit by Eq. (g) agrees very well in the whole range 
explored except for the region of the bump ascribed to 
the BP, which is not accounted for in the fitting func- 
tion. Since we can carry out a precise line shape analysis 
by using only two relaxation processes, where the long 
time part is perfectly described by the a relaxation of 
the MCT, we can test the predictions of the MCT which 
concern the values for the Kohlrausch exponent f3 and 
the long relaxation time 77 . 
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FIG. 11. The free water contribution to the ISF for 
Q = 2.25 A -1 along the z direction at different tempera- 
tures, from T = 298 K (bottom curve) to T = 210 K (top 
curve). The curves are fitted by Eq. (H) (long dashed lines). 
It is evident that the bump, the region of which is excluded 
from the fit, increases with decreasing temperature. 
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FIG. 12. Fit parameters. Top: log-log plots of the inverse 
relaxation time t^ 1 as a function of T — Tc along the z 
direction (left) and the xy direction (right). The full lines 
are power law fits given by 1/t z ~ (T — 167. 6) 3 ' 8 (left) and 
1/Txy ~ (T- 170.4) 3 ' 4 (right). Bottom: T dependence of the 
Kohlrausch exponents /3 along the z (left) and xy directions 
(right). 

The values of T\ and /3 obtained from the fits of the ISF 
at different temperatures are shown in Fig. |l2| for both xy 
and z directions. The validity of the asymptotic scaling 
law predicted by MCT implies that .the stretching expo- 
nent (3 is temperature independent .O For T > Tc one 
expects f3(T — ► Tc) ~ A) < 1- In our case the exponent 
[3 approaches an asymptotic value around 0.5 in direc- 
tion and 0.6 along z as shown in the lower part of Fig. [l^. 
A similar value is found for bulk water. In the panels in 
the upper part we show that the relaxation time satisfies 
the MCT prediction Eq. (^) . In the z direction (panel on 
the left) we estimate a temperature of structural arrest 
Tc ~ 167.6 K with an exponent 7 w 3.8. Along the xy 
direction (panel on the right) the Tc is slightly higher 
T c w 170.4 K with 7 w 3.4. 

This is in agreement with the idea that an enhance- 
ment of the free volume speeds up the dynamics leading 
to a decrease of the glass transition temperature relative 
to the bulk. At ambient pressure SPC/E bulk water un- 
dergoes a_kinetic glass transition at T = 186.3 K with 
7 = 2.29.N 



B. The ISF at different Q 

From the fits of the ISF at different wavelengths Q, 
shown for example along the z direction in Fig. [l3| for 
T = 240 K, we obtain values of and 77 as a function 
of Q. Such values for the different temperatures are re- 
ported in Fig. [l4| for the z direction and in Fig. |l5| for the 
xy direction. The behavior of the parameters is compati- 
ble with the MCT. As shown in the insets of both figures 
the Kohlrausch exponent (3 starts from values close to 
1 at low Q and high temperatures and reaches a com- 
mon plateau at large Q. This behavior is in agreement 
with the interpretation of the kinetic glass transition by 
MCT in terms of the cage effect, since for Q^ 1 larger 
than the cage size we expect to find that the stretching 
of the exponential at long time Eq. (0) is less relevant. 
Thus the Kohlrausch exponent is expected to go to 1 
for large Q^ 1 when the system enters the normal Debye 
regime of diffusion in the long time decay region of the 
ISF. Of course this limit is reached at higher values of Q 
at higher temperatures, where the Debye regime is more 
easily established, as seen in the insets of Fig. 
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and If 



The oscillation in the plateau region does not allow for 
the calculation of the MCT parameters a and b related 
to the power laws of approach to and departure from the 
plateau. It is however shown that for large Q values the 
Kohlrausch exponent coincides with the scaling exponent 
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bJZ3 We therefore obtain b w 0.4 in this system for both 
the xy and the z directions. 
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FIG. 13. ISF along the z direction of the free water for 
half hydration and T — 240 K for several values of Q from 
Q = 1.8 A" 1 (top curve) to Q = 3 A -1 (bottom curve). All 
curves are fitted by Eq. (@) (long dashed lines), excluding the 
region around the bump. 




FIG. 14. Relaxation time n and exponent f3 (in the in- 
set) of free water as extracted from the fits of the long time 
relaxation behavior of the ISF at different Q along the xy 
direction (Eq. (^)). The continuous lines with slope 1 have 
been inserted as guidelines. Different symbols represent differ- 
ent temperatures; temperature decreases from top to bottom. 
Open symbols in the inset are for the same temperature as 
the corresponding full symbols in the figure. 

rf 1 shows a linear behavior as function of Q 2 in Fig. [l4| 
and [l5[ where the continuous lines which have slope 1 are 
plotted as guidelines. This behavior is common for glass 
formers close to the kinetic glass transition and has been 



found for instance in glycerol. In the interpretation we 
must take into account that r; measures the timescale of 
the final relaxation of the cage with the subsequent decay 
of the self correlation function of the densities over the 
length scale Q~ x and the diffusion of the molecule over 
distances of the order Q~ x . In the diffusive regime valid 
for large Q~ x we expect t, ~ Q 2 . At variance with 
bulk water we do not find any crossover of the behaviau 
of rf from a Q 2 to a Q proportionality at high Q,ca 
except, perhaps, at the lowest temperature. 

In Fig. there is a clear evidence of a bump at 0.7 ps 
which ispreserved at the different temperatures as seen 
in Fig. [ll] for Q max — 2.25 A -1 . The position of the 
bump is independent of Q, which is clearly seen in Fig. [l3| 
at T = 240 K and which is also found at all other ex- 
plored temperatures. As mentioned in Sec. VI A, similar 
effects have been nbsftrvfid, in computer simulations of 
other glass formersEScJCao and have been attributed to 
the existence of a Boson peak (BP). The BP is present in 
many glasses as a consequence of an excess of vibrational 
modes and can be detected even ia computer simulation 
in spite of the finite size effects. HEj We do not want to 
enter here in a detailed discussion of this behavior, but 
we recall that the BP is considered as a precursor of a 
glass transition when it appears in the liquid state. We 
notice that as seen in Fig. the BP feature appears 
clearly for the inner shells while it is much less evident in 
the outer layer contribution. While in some respect the 
substrate seems to enhance the effect relative to bulk wa- 
tcicJ, it appears that water in the outer shells is not able 
to sustain the vibrations which are considered responsi- 
ble of the overshoot, possibly because we keep the atoms 
of the silica substrate rigid during the simulation. We 
intend to discuss this point more extensively in a future 
paper. 




FIG. 15. Like Fig. |lj, but along the z direction. 
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VII. COMPARISON WITH EXPERIMENTAL 
DATA 

The technique that allows for the best comparison 
with our MD results is quasielastic neutron scattering 
(QENS). In fact the self dynamic structure factor as mea- 
sured by incoherent QENS is the time Fourier transform 
of the ISF. The first data appeared in literature on water- 
in- Vycor were well fitted with a jump-diffusion model in 
a confined medium.O This model was giving reasonable 
relaxation times, but the size of the confining region was 
unrealistic, 5 A. An analysis-based on MCT of those data 
was successively proposed. The data were also well 
fitted and the results were consistent with the appear- 
ance upon supercooling of a slow dynamics due to the 
approach to the kinetic glass transition. The Q depen- 
dence of the inverse of the structural relaxation time ap- 
peared for the 50% hydration level upon supercooling, 
as a .power law with an exponent ranging from 2.51 to 
2.06.EI The fits were however performed on a limited Q 
range, approximately from 0.4 to 1.3 A -1 . These find- 
ings have been confirmed by a recent QENS experiment 
where a full set of high resolution data on water confined 
in Vycor glass has been analyzed for different hydration 
levels upon supercoolingEd The experimental exponents 
obtained from the fits to a limited range of Q appear to be 
slightly different from ours. In fact we find, see Figs. |l4|- 
|l5| , a value which is approximately 2 at high temperatures 
and slightly decreases for the lowest temperature. How- 
ever an analysis performed on a more extended Q value 
data, up to 1.8 A" 1 showed an exponent which is 2 at 
high temperature and decreases at lower temperatures, 
similar to what we found in the present simulation.^ 

It is important to notice that the QENS response func- 
tion is able to probe only the mobile water, namely the 
water above Tc, because of experimental resolution lim- 
itations. Therefore the experimental data appear so far 
compatible with our findings. 



VIII. SUMMARY AND CONCLUSIONS 

In this paper we presented the results obtained by 
Molecular Dynamics on the single particle dynamics of 
SPC/E water confined in a silica pore, which represents 
an approximation to the system of water confined in Vy- 
cor. This system can also be considered as representative 
of the more general case of water confined in a hydrophilic 
nanopore. 

Through the study of dynamic time correlators such as 
the mean square displacement (MSD) and the self part of 
the intermediate scattering factor (ISF), we investigated 
the onset of slow dynamics when the hydration level in 
the pore is decreased or when the liquid is supercooled. 

At ambient temperature we found that upon decreas- 
ing the hydration level a flattening of the MSD takes 



place after the initial ballistic regime and before the on- 
set of the stochastic diffusion. In correspondence a two 
step relaxation behavior appears in the ISF, which dis- 
plays a strongly non-exponential decay at long time. In 
this way the single particle dynamics has the signatures 
of the presence of a cage effect as predicted by Mode 
Coupling Theory (MCT). 

Apparently the effect of hydration is similar to the ef- 
fect of supercooling in bulk water, but the long time re- 
laxation of the ISF cannot be fitted by the stretched ex- 
ponential predicted by the MCT. By performing a layer 
analysis it becomes possible to separate the contribution 
to the dynamics of the water molecules close to the sub- 
strate (bound water) from the contribution coming from 
the molecules which are in the internal layers (free wa- 
ter). The subset of molecules belonging to the bound 
water shows density oscillations and appears to be in a 
glassy state with very low mobility already at ambient 
conditions. The free water shows a dynamical behavior 
typical of glass forming liquids. This behavior can be 
accounted for by the idealized MCT in analogy withjjijfc. 
vious findings on SPC/E water in the bulk phase. L3C3'cZl 
We would like to stress, however, that in the present 
work we have considered a much more complex system 
than the bulk, where an additional anisotropic effect and 
a perturbation by the substrate is present. 

The fact that MCT could be used in this framework 
as a unifying theoretical approach is highly relevant as a 
guideline for the systematic study of the important phe- 
nomenology of confined and interfacial water. 

We show evidence that the hydration level together 
with the hydrophilicity of the substrate play, through 
the layering effect, an important role in determining the 
dynamical properties of confined water. As far as the 
anisotropy of the system is concerned, we did not ob- 
serve large differences in the behavior along the confined 
direction (xy) and along the pore axis (z), where periodic 
boundary conditions are applied, apart from the fact that 
the mobility along the xy direction is smaller than in the 
z direction. 

By supercooling confined water at half hydration we 
have found that the subset of free water undergoes a ki- 
netic glass transition, which can be described in terms of 
the MCT. The layer analysis makes possible the calcula- 
tion of various important quantities related to the glass 
transition, like the temperature of structural arrest, Tc, 
and the critical exponent 7. Similar results have been ob- 
tained, from a preliminary analysis of the full hydration 
case.o 

MCT has been largely used in the interpretation of the 
phenomenology of liquids in the supercooled phase. In 
the recent attempts to develop a unified first-principles 
theory for the glass transition it has been recognized that 
MCT keeps its validity in describing the liquid just above 
the temperature To It is time to test this theory for 
confined fluids, where, in spite of the experimental and 
computer simulation efforts, the phenomenology of the 
dynamics of the confined and interfacial liquids upon su- 
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percooling is still unclear. We have shown evidence that 
the MCT can be applied to the interpretation of the dy- 
namical behavior of supercooled confined water, if the 
layering effects determined by the interaction with the 
hydrophilic substrate are carefully taken into account. 
We believe that the results obtained in this work are of 
great importance in understanding the phenomenology 
of the mobility of fluids in nanopores. 
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